/* MolObj.m - Copyright 1993  Steve Ludtke */
/* This is the central object in the program. Almost everything */
/* goes through this object on some level */

#import "Mtypes.h"
#import "MolObj.h"
#import "Inspector.h"
#import "D3View.h"
#import "AminoView.h"

#import <stdio.h>
#import <string.h>
#import <libc.h>
#import <math.h>
#import <appkit/appkit.h>
#import <dpsclient/psops.h>

#define NSTEP 50
#define NSTEP0 200
#define swapa(n1,n2) { tmp=mol[0].am[n1]; mol[0].am[n1]=mol[0].am[n2]; mol[0].am[n2]=tmp; }

/* vector operations to simplify math, boy, c++ would be nice ... */
float dot(float *a,float *b);
void cross(float *a,float *b,float *c);
void sub(float *a,float *b,float *c);
float len(float *a);
void count (Molecule *mol,int *cnt, int n1,short lev,short maxlev);

void tagatoms(Atom *atom,int n);
float frand(float lo, float hi);
extern ACID acid[AADEF];

@implementation MolObj

-init
{
FILE *in;
char *bpath;
char s[80];
int i;

[super init];

srandom(time(0));

/* find where the app is located */
bpath=(char *)[[NXBundle mainBundle] directory];
/* read the covalent.rad file with atom preferences */
sprintf(s,"%s/Library/MolViewer.dat",getenv("HOME"));
in=fopen(s,"r");
if (in==NULL) {
	sprintf(s,"%s/MolViewer.dat",bpath);
	in=fopen(s,"r");
	if (in==NULL) { printf("No .rad\n"); exit(0); }
}
for (i=0; fgets(s,79,in)!=NULL&&i<NEL; i++) {
	elinfo[i].color[0]=elinfo[i].color[1]=elinfo[i].color[2]=.7;
	elinfo[i].name[0]=s[0];
	if (s[1]!='\t') elinfo[i].name[1]=s[1];
	else elinfo[i].name[1]=' ';
	sscanf(&s[2],"%f %f %f %f",&elinfo[i].arad,&elinfo[i].color[0],&elinfo[i].
		color[1],&elinfo[i].color[2]);
/* the 1.3 is user adjustable after the program has been started */
	elinfo[i].rad=1.3*elinfo[i].arad;
	elinfo[i].image=nil;
}
fclose(in);

/* initialize arrays */
for (i=0; i<MAXMOL; i++) { mol[i].verts=NULL; mol[i].nverts=NULL; }

/* load shapes for quick spacefilling mode */
sprintf(s,"%s/h.tiff",bpath);
elinfo[0].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
sprintf(s,"%s/c.tiff",bpath);
elinfo[5].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
sprintf(s,"%s/n.tiff",bpath);
elinfo[6].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
sprintf(s,"%s/o.tiff",bpath);
elinfo[7].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
sprintf(s,"%s/p.tiff",bpath);
elinfo[14].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];
sprintf(s,"%s/s.tiff",bpath);
elinfo[15].image=[[[NXImage alloc] initFromFile:s] setScalable:YES];

nmol=0;
stepmode=1;
cstep=1.0;
return self;
}

/* we're the App's delegate, so tell everyone else that we're ready to go */
-appDidInit:sender
{
[molView appDidInit:sender];
[molView setData:mol :nmol :elinfo];
[inspector startup:mol :nmol :elinfo];
return self;
}

/* read a MolViewer file (format in transition, don't use) */
-readMol:sender
{
id panel;
FILE *in;
int i,j,k;
float x,y,z,x1,y1,z1;
char c[2],s[101];

x=y=z=0.0;

panel=[[OpenPanel new] allowMultipleFiles:NO];
if (![panel runModalForTypes:NULL]) return self;
in=fopen([panel filename],"r");
if (in==NULL) return self;

[self remMol:self];
while (fgets(s,100,in)!=NULL);
if (sscanf(s," %d",&mol[nmol].numa)!=1 || mol[nmol].numa>5000)
	{ fclose(in); return self; }
mol[nmol].numa++;
mol[nmol].atom=malloc(sizeof(Atom)*(mol[nmol].numa+1));
mol[nmol].coord=malloc(sizeof(RtPoint)*(mol[nmol].numa+1));
rewind(in);

for (i=0; fgets(s,100,in)!=NULL; i++) {
	if (s[0]=='#') { i=0; continue; }
	sscanf(s,"%d %c%c %f %f %f",&j,c,&c[1],&x1,&y1,&z1);
	if (c[1]=='\t') c[1]=' ';

	mol[nmol].coord[j][0]=x1;
	mol[nmol].coord[j][1]=y1;
	mol[nmol].coord[j][2]=z1;
	mol[nmol].atom[j].numb=0;
	mol[nmol].atom[j].cnum=i;	/* atom # within residue (for proteins) */
	mol[nmol].atom[j].tag=mol[nmol].atom[j].sel=0;

	for (k=0; k<NEL; k++) if (c[0]==elinfo[k].name[0]&&c[1]==elinfo[k].name[1]) break;
	if (k==NEL) printf("Can't find atom:%c\n",c[0]);
	mol[nmol].atom[j].anum=k;

	x+=x1;	/* keep running total for finding average x,y,z */
	y+=y1;
	z+=z1;
}
fclose(in);
j++;
x/=(float)j;
y/=(float)j;
z/=(float)j;

/* move origin to center of molecule */
for (i=0; i<mol[nmol].numa; i++) {
	mol[nmol].coord[i][0]-=x;
	mol[nmol].coord[i][1]-=y;
	mol[nmol].coord[i][2]-=z;
}

mol[nmol].numlb=0;
strcpy(s,[panel filename]);
s[strlen(s)-3]='b';
s[strlen(s)-2]='n';
s[strlen(s)-1]='d';
in=fopen(s,"r");
if (in==NULL) { printf("Can't open bnd\n"); nmol++; 
	[molView setData:mol :nmol :elinfo]; return self; }

while (fgets(s,100,in)!=NULL) {
	if (s[0]=='L') mol[nmol].numlb++;
}
mol[nmol].lb=malloc(sizeof(struct LBOND)*mol[nmol].numlb);
rewind(in);

i=j=0;
while (fgets(s,100,in)!=NULL) {
	if (s[0]=='L') {
		sscanf(&s[1],"%d %d %f %f",&mol[nmol].lb[i].n1,&mol[nmol].lb[i].n2, 
			&mol[nmol].lb[i].d,&mol[nmol].lb[i].w);
		if (mol[nmol].lb[i].n1>mol[nmol].numa-1 || 
			mol[nmol].lb[i].n2>mol[nmol].numa-1) i--;
		i++;
	}
}
[self asnBnd:nmol];
nmol++;
[inspector newMol];
[self update:U_TOTAL];
return self;
}

/* read a Protein Data Bank molecule */
-readPdb:sender
{
id panel;
FILE *in;
int i,j,k,nb,b[8];
float x,y,z;
char s[101];
RtPoint *a1,*a2;

panel=[[OpenPanel new] allowMultipleFiles:NO];
if (![panel runModalForTypes:NULL]) return self;
in=fopen([panel filename],"r");
if (in==NULL) return self;

/* get rid of the current molecule */
[self remMol:self];
i=nb=0;
k=1;
x=y=z=0;
/* first pass, count atoms and bonds */
while (fgets(s,100,in)!=NULL) {
/* still uses CONECT records to generate bonds. Needs to be replaced */
	if (strncmp(s,"CONECT",6)==0) {
		s[70]=0;
		j=sscanf(&s[7]," %d %d %d %d %d %d %d %d",&k,&k,&k,&k,&k,&k,&k,&k);
		nb+=j-1;
	}
	if (k && (strncmp(s,"ATOM",4)==0 || strncmp(s,"HETATM",6)==0)) {
		sscanf(&s[7]," %d",&i);
	}
/* only reads the first structure in the file, TER terminates reading */
	if (strncmp(s,"TER",3)==0) k=0;
}
rewind(in);
mol[nmol].numa=i;
mol[nmol].atom=malloc(sizeof(Atom)*(i+1));
mol[nmol].coord=malloc(sizeof(RtPoint)*(i+1));
mol[nmol].numlb=nb;
if (nb!=0) mol[nmol].lb=malloc(sizeof(struct LBOND)*(nb+1));

i=nb=0;
k=0;
/* ok, this time actually read the data */
while (fgets(s,100,in)!=NULL) {
	if (strncmp(s,"CONECT",6)==0) {
		s[70]=0;
		j=sscanf(&s[7]," %d %d %d %d %d %d %d %d",b,b+1,b+2,b+3,b+4, 
			b+5,b+6,b+7)-1;
		for (i=0; i<j; i++) {
			a1=(RtPoint *)&mol[nmol].coord[b[i]-1][0];
			a2=(RtPoint *)&mol[nmol].coord[b[i+1]-1][0];
			if (b[i]>mol[nmol].numa || b[i+1]>mol[nmol].numa) continue;
			mol[nmol].lb[nb].n1=b[i]-1;
			mol[nmol].lb[nb].n2=b[i+1]-1;
			mol[nmol].lb[nb].w=1.0;
			mol[nmol].lb[nb].d=sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+ 
				SQR(a1[2]-a2[2]));
			nb++;
		}
	continue;
	}
	if (k) continue;
	if (strncmp(s,"TER",3)==0) k=1;
	if (strncmp(s,"ATOM",4)!=0 && strncmp(s,"HETATM",6)!=0) continue;
	sscanf(&s[6]," %d",&i);
	i--;
/* try to find this atom in the atomic table */
	for (j=0; j<NEL; j++) 
		if(s[12]==elinfo[j].name[0]&&s[13]==elinfo[j].name[1]) break;
	if (j==NEL) {
		for (j=0; j<NEL; j++) if(s[13]==elinfo[j].name[0]) break;
	}
	if (j==NEL) printf("Can't find atom:%c%c\n",s[12],s[13]);
	mol[nmol].atom[i].anum=j;
	mol[nmol].atom[i].numb=0;
	mol[nmol].atom[i].tag=mol[nmol].atom[i].sel=0;
	sscanf(&s[30]," %f %f %f",&mol[nmol].coord[i][0], 
		&mol[nmol].coord[i][1],&mol[nmol].coord[i][2]);
/* find the "average" position for centering */
	x+=mol[nmol].coord[i][0];
	y+=mol[nmol].coord[i][1];
	z+=mol[nmol].coord[i][2];
}
mol[nmol].numlb=nb;

x/=(float)mol[nmol].numa;
y/=(float)mol[nmol].numa;
z/=(float)mol[nmol].numa;

/* move origin to center of molecule */
for (i=0; i<mol[nmol].numa; i++) {
	mol[nmol].coord[i][0]-=x;
	mol[nmol].coord[i][1]-=y;
	mol[nmol].coord[i][2]-=z;
}

[self asnBnd:nmol];
nmol++;
[inspector newMol];
[self update:U_TOTAL];
return self;
}

/* read an Alchemy format file */
-readAlch:sender
{
id panel;
FILE *in;
int i,j,k,l;
float x,y,z,x1,y1,z1;
float *a1,*a2;
char c[4],s[101];

x=y=z=0.0;

panel=[[OpenPanel new] allowMultipleFiles:NO];
if (![panel runModalForTypes:NULL]) return self;
in=fopen([panel filename],"r");
if (in==NULL) return self;

[self remMol:self];
/* Ah Ha! Alchemy has a header with an atom/bond count */
fgets(s,100,in);
sscanf(s," %d",&mol[nmol].numa);
sscanf(&s[13]," %d",&mol[nmol].numlb);
mol[nmol].atom=malloc(sizeof(Atom)*(mol[nmol].numa+1));
mol[nmol].coord=malloc(sizeof(RtPoint)*(mol[nmol].numa+1));

for (i=0; i<mol[nmol].numa; i++) {
	fgets(s,100,in);
	if (s[0]=='#') { i=0; continue; }
	sscanf(s," %d %c%c%c%c %f %f %f",&j,c,&c[1],&c[2],&c[3],&x1,&y1,&z1);
	j--;

	mol[nmol].coord[j][0]=x1;
	mol[nmol].coord[j][1]=y1;
	mol[nmol].coord[j][2]=z1;
	mol[nmol].atom[j].numb=0;
	mol[nmol].atom[j].cnum=i;	/* atom # within residue (for proteins) */

/* try to identify the atoms. Free bonds become hydrogens */
	if (c[0]=='~') c[0]='H';
	for (k=0; k<NEL; k++) {
		if ((c[1]>='0'&&c[1]<='9')||c[1]==' '||c[2]!=' ') {
			if (c[0]==elinfo[k].name[0]&&elinfo[k].name[1]==' ') break;
		}
		else if (c[0]==elinfo[k].name[0]&&c[1]==elinfo[k].name[1]) break;
	}
	if (k==NEL) printf("Can't find atom:%c\n",c[0]);
	mol[nmol].atom[j].anum=k;
	mol[nmol].atom[j].type[0]=c[1];
	mol[nmol].atom[j].type[1]=c[2];
	mol[nmol].atom[j].type[2]=c[3];
	mol[nmol].atom[j].tag=mol[nmol].atom[j].sel=0;
	x+=x1;	/* keep running total for finding average x,y,z */
	y+=y1;
	z+=z1;
}
j++;
x/=(float)j;
y/=(float)j;
z/=(float)j;

/* move origin to center of molecule */
for (i=0; i<mol[nmol].numa; i++) {
	mol[nmol].coord[i][0]-=x;
	mol[nmol].coord[i][1]-=y;
	mol[nmol].coord[i][2]-=z;
}

mol[nmol].lb=malloc(sizeof(struct LBOND)*mol[nmol].numlb);

/* now read the bonds */
j=0;
for (i=0; i<mol[nmol].numlb; i++) {
	fgets(s,100,in);
	sscanf(s,"%d %d %d",&j,&k,&l);
	mol[nmol].lb[i].n1=k-1;
	mol[nmol].lb[i].n2=l-1;
	a1=(float *)&mol[nmol].coord[k-1][0];
	a2=(float *)&mol[nmol].coord[l-1][0];
	mol[nmol].lb[i].w=1.0;
	mol[nmol].lb[i].type=s[19];
	x=sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+SQR(a1[2]-a2[2]));
	if (x>2.2) x=1.8;
	mol[nmol].lb[i].d=x;
}

/*for (i=0; i<mol[0].numa; i++) {
	mol[0].coord[i][1]*=.5;
	mol[0].coord[i][0]=frand(-4.0,4.0);
	mol[0].coord[i][2]=frand(-4.0,4.0);
}*/

[self asnBnd:nmol];
nmol++;
fclose(in);
[inspector newMol];
[self update:U_TOTAL];
return self;
}

/* delete the current molecule from memory */
-remMol:sender
{
if (nmol==0) return self;
nmol--;
if (mol[nmol].numa!=0) free(mol[nmol].atom);
if (mol[nmol].numlb!=0) free(mol[nmol].lb);
if (mol[nmol].verts!=NULL) { free(mol[nmol].verts); mol[nmol].verts=NULL; }
if (mol[nmol].nverts!=NULL) { free(mol[nmol].nverts); mol[nmol].nverts=NULL; }
mol[nmol].numa=mol[nmol].numlb=0;
[self update:U_TOTAL];
return self;
}

/* something has changed, so make sure all the displays are correct */
-update:(int)level
{
if (level>=U_BONDS) [self setupMol];
[molView setData:mol :nmol :elinfo];
[inspector update:nmol :level];
return self;
}

/* generates a data table to speed up drawing */
-setupMol
{
int i,j,k;

for (i=0; i<nmol; i++) {
	k=mol[i].numlb;
	if (mol[i].verts!=NULL) free(mol[i].verts);
	if (mol[i].nverts!=NULL) free(mol[i].nverts);
	mol[i].verts=malloc(sizeof(RtInt)*2*(k+1));
	mol[i].nverts=malloc(sizeof(RtInt)*(k+1));
	for (j=0; j<k; j++) {
		mol[i].nverts[j]=2;
		mol[i].verts[j*2]=mol[i].lb[j].n1;
		mol[i].verts[j*2+1]=mol[i].lb[j].n2;
	}
}
return self;
}

/* Write an Alchemy format file */
-dump:sender
{
int i;
static id savePanel=nil;
FILE *out;
char s[80];

if (savePanel==nil) savePanel=[SavePanel new];
[savePanel setRequiredFileType:"mol"];
  
if (![savePanel runModal]) return self;
out=fopen([savePanel filename],"w");

fprintf(out,"%5d ATOMS, %5d BONDS,     0 CHARGES, UNK\n", 
			mol[0].numa,mol[0].numlb);
for (i=0; i<mol[0].numa; i++) {
	fprintf(out,"%5d %c%c%c%c  %8.4f %8.4f %8.4f     0.0000\n", i+1,elinfo[mol[0].atom[i].anum].name[0],mol[0].atom[i].type[0],mol[0].atom[i].type[1],mol[0].atom[i].type[2],mol[0].coord[i][0],mol[0].coord[i][1], mol[0].coord[i][2]);
}

for (i=0; i<mol[0].numlb; i++) {
	if (mol[0].lb[i].type=='S') strcpy(s,"SINGLE");
	else strcpy(s,"DOUBLE");
	fprintf(out,"%5d %5d %5d  %s\n",i+1,mol[0].lb[i].n1+1,mol[0].lb[i].n2+1,s);
}
fclose(out);
return self;
}

/* copy bond table to redundant bond information in Atom structure */
-asnBnd:(int)nm
{
int i,j;
int lb,n1,n2;
float l;

lb=mol[nm].numlb;

for (i=0; i<lb; i++) {
	n1=mol[nm].lb[i].n1;
	n2=mol[nm].lb[i].n2;
	l=mol[nm].lb[i].d;

	j=mol[nm].atom[n1].numb++;
	if (j>3) { j=3; /*printf("Too many bonds\n");*/ }
	mol[nm].atom[n1].bnd[j]=n2;
	mol[nm].atom[n1].bndl[j]=l;

	j=mol[nm].atom[n2].numb++;
	if (j>3) { j=3; /*printf("Too many bonds\n");*/ }
	mol[nm].atom[n2].bnd[j]=n1;
	mol[nm].atom[n2].bndl[j]=l;
}

for (i=0; i<mol[nm].numa; i++) {
	for (j=mol[nm].atom[i].numb; j<4; j++) mol[nm].atom[i].bnd[j]=-1;
}

if (nm==0) chi=[self calcChi:self];
return self;
}

/* locate and identify amino acids in the current molecule */
-findAA:sender
{
int i,j,k,l,h=0,n,na;
int cnt[16];
struct AMINO tmp;

/* This is difficult to follow, but you'll just have to wade through it */
/* if you need to change something. */
na=0;
for (i=0; i<mol[0].numa; i++) {
	if (mol[0].atom[i].anum!=6 || mol[0].atom[i].type[0]!='A') continue;
	mol[0].am[na].n=i;			/* assign N */
	mol[0].atom[i].tag=1;
	mol[0].am[na].cn= -1;
								/* assign h,c1, and cn */
	for (j=0; j<3; j++) {
		k=mol[0].atom[i].bnd[j];
/* h is attached to n */	
		if (mol[0].atom[k].anum==0) mol[0].am[na].h=k;
/* cn has an O attached, c1 has an H and the sidechain */
		if (mol[0].atom[k].anum==5) {
			if (mol[0].atom[k].type[0]=='3') mol[0].am[na].c1=k;
			else if (mol[0].atom[k].type[0]=='2') mol[0].am[na].cn=k;
			else printf("Bad carbon on Nam\n");
		}
	}
								/* assign cn2 */
	k=mol[0].am[na].cn;
	if (k==-1) mol[0].am[na].cn2= -1;
	else {
/* cn2 is attached to cn and an o */
		for (j=0; j<mol[0].atom[k].numb; j++) {
			n=mol[0].atom[k].bnd[j];
			if (mol[0].atom[n].anum==5 && mol[0].atom[n].type[0]=='3') 
				mol[0].am[na].cn2=n;
		}
	}
								/* assign c2 */
	k=mol[0].am[na].c1;
	mol[0].atom[k].tag=1;
/* c2 is attached to c1, h, and a sidechain */
	for (j=0; j<mol[0].atom[k].numb; j++) {
		n=mol[0].atom[k].bnd[j];
		if (mol[0].atom[n].anum==5 && mol[0].atom[n].type[0]=='2') 
			mol[0].am[na].c2=n;
	}
								/* assign o */
	k=mol[0].am[na].c2;
	for (j=0; j<mol[0].atom[k].numb; j++) {
		n=mol[0].atom[k].bnd[j];
		if (mol[0].atom[n].anum==7) mol[0].am[na].o=n;
	}
								/* assign sc */
	k=mol[0].am[na].c1;
	for (j=0; j<4; j++) {
		n=mol[0].atom[k].bnd[j];
		if (n==mol[0].am[na].c2||n==mol[0].am[na].cn) continue;
		if (mol[0].atom[n].anum==0) { h=n; continue; }
		if (mol[0].atom[n].anum==5) { mol[0].am[na].sc=n; break; }
	}
/* if glycine, an h is chosen randomly */	
	if (j==4) mol[0].am[na].sc=h;

	na++;
	if (na==MAXAA) { na--; printf("Max # of amino acids exceeded\n"); }
}	

mol[0].numaa=na;

for (i=0; i<mol[0].numaa; i++) {
/* calculate the dihedrals */
	mol[0].am[i].sel=0;
	mol[0].am[i].omega=[self dihedral:mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].cn :mol[0].am[i].cn2];
	mol[0].am[i].psi=[self dihedral:mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].c2 :mol[0].am[i].o]-M_PI;
	mol[0].am[i].phi=[self dihedral:mol[0].am[i].c2 :mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].h]-M_PI;
	mol[0].am[i].anum=0;
	if (mol[0].am[i].psi< -M_PI) mol[0].am[i].psi+=M_PI;
	if (mol[0].am[i].phi< -M_PI) mol[0].am[i].phi+=M_PI;

/* try to identify the individual amino acids by counting and a little */
/* chain testing when the result is ambigous. Exotic amino acids */
/* might be incorrectly assigned. Proline is also incorrectly assigned. */ 
	for (j=0; j<16; j++) cnt[j]=0;
	j=mol[0].am[i].sc;
	if (mol[0].atom[j].anum==0) mol[0].am[i].anum=9;
	else {
		count(mol,cnt,mol[0].am[i].sc,0,15);	
		for (k=0; k<AADEF; k++) {
			if (acid[k].nc==cnt[5]&&acid[k].ns==cnt[15]&&acid[k].no==cnt[7]&& acid[k].nn==cnt[6]) break;
		}
		if (k==AADEF) { mol[0].am[i].anum=0;
			printf("C=%d S=%d O=%d N=%d\n",cnt[5],cnt[15],cnt[7],cnt[6]);
		}
		else mol[0].am[i].anum=k;
	}
	if 	(mol[0].am[i].anum==3) {
		l=mol[0].am[i].sc;
		cnt[5]=0;
		for (k=0; k<mol[0].atom[l].numb; k++) if (mol[0].atom[mol[0].atom[l].bnd[k]].anum==5) cnt[5]++;
		if (cnt[5]==3) mol[0].am[i].anum=4;
	}
}
for (i=0; i<mol[0].numa; i++) mol[0].atom[i].tag=0;

/* sort amino acids so N terminus is FIRST */
for (i=1; i<mol[0].numaa; i++) if (mol[0].am[i].cn==-1) swapa(0,i);
for (i=1; i<mol[0].numaa; i++) {
	for (j=i+2; j<mol[0].numaa; j++) {
		if (mol[0].am[i-1].c2==mol[0].am[j].cn) swapa(i,j);
	}
}
return self;
}

/* recursive routine to decend the sidechain and count the number of */
/* atoms of each type */
void count (Molecule *mol,int *cnt, int n1,short lev,short maxlev)
{
int i,j;

if (lev>maxlev||mol[0].atom[n1].tag) return;
j=mol[0].atom[n1].anum;
if (j<16) cnt[j]++;
mol[0].atom[n1].tag=1;
for (i=0; i<mol[0].atom[n1].numb; i++) {
	count(mol,cnt,mol[0].atom[n1].bnd[i],lev+1,maxlev);
}
return;
}

/* calculates a "chi" for fitting (not used yet ...) */
-(float)calcChi:sender
{
int i;
int lb;
float l,ch;
float *a1,*a2;

lb=mol[0].numlb;
ch=0.0;
for (i=0; i<lb; i++) {
	a1=(float *)&mol[0].coord[mol[0].lb[i].n1];
	a2=(float *)&mol[0].coord[mol[0].lb[i].n2];
	l=mol[0].lb[i].d;
	ch+=fabs(l-sqrt(SQR(a1[0]-a2[0])+SQR(a1[1]-a2[1])+SQR(a1[2]-a2[2])));
}
ch/=lb;
return ch;
}

/* fitting stuff ... */
-minStep:sender
{
int i,j,k,m,c,c2;
float x,y,z,d,d2,v1,v2;
float x1,y1,z1;
float *a0,*a1;
short *ba;

c2=0;
if (stepmode==0) {
	for (i=0; i<NSTEP0; i++) {
		x1=y1=z1=0.0;
		j=random()%mol[0].numa;
		k=mol[0].atom[j].numb;
		a0=(float *)&mol[0].coord[j];
		for (m=0; m<k; m++) {
			a1=(float *)&mol[0].coord[mol[0].atom[j].bnd[m]];
			x=a1[0]-a0[0];
			y=a1[1]-a0[1];
			z=a1[2]-a0[2];
			d2=sqrt(SQR(x)+SQR(y)+SQR(z));
			d=(1.0/mol[0].atom[j].bndl[m]-1.0/d2);
			x1+=x*d;
			y1+=y*d;
			z1+=z*d;
		}
		x1/=(float)k;
		y1/=(float)k;
		z1/=(float)k;
		a0[0]+=x1;
		a0[1]+=y1;
		a0[2]+=z1;
	}
	chi=[self calcChi:self];
	c2=100;
}
else if (stepmode==1) {
	for (i=0; i<NSTEP; i++) {
		j=random()%mol[0].numa;
		m=1;
		c=0;

		ba=mol[0].atom[j].bnd;
		a0=(float *)&mol[0].coord[j];
		v1=0.0;
		for (k=0; k<mol[0].numa; k++) {
			if (k==j||k==ba[0]||k==ba[1]||k==ba[2]||k==ba[3]) continue;
			a1=(float *)&mol[0].coord[k];
			d=sqrt(SQR(a0[0]-a1[0])+SQR(a0[1]-a1[1])+SQR(a0[2]-a1[2]));
			v1+=1.0/d;
		}

		while (m&&c++<30) {
			x1=frand(-cstep,cstep);
			y1=frand(-cstep,cstep);
			z1=frand(-cstep,cstep);

			v2=0.0;
			a0[0]+=x1;
			a0[1]+=y1;
			a0[2]+=z1;

			m=0;			
			for (k=0; k<mol[0].numa; k++) {
				if (k==j||k==ba[0]||k==ba[1]||k==ba[2]||k==ba[3]) continue;
				a1=(float *)&mol[0].coord[k];
				d=sqrt(SQR(a0[0]-a1[0])+SQR(a0[1]-a1[1])+SQR(a0[2]-a1[2]));
				v2+=1.0/d;
			}
/*			if (v1<v2) m=1; else m=0;
			if (m) { a0[0]-=x1; a0[1]-=y1; a0[2]-=z1; continue; } */
			d=v1/(mol[0].numa*6.0);
			v2/=mol[0].numa*6.0;
			x=[self calcChi:self];
			if (x+v2>chi+d) { m=1; a0[0]-=x1; a0[1]-=y1; a0[2]-=z1; continue; }
			chi=x;
		}
		if (c<=30) c2++;
	}
	[chiDevDis setFloatValue:chi];
}
[c2Dis setIntValue:c2];
if (c2<3) {
	cstep*=.9;
	[csDis setFloatValue:cstep];
}
return self;
}

/* more fitting stuff */
-resetCs:sender
{
cstep=1.0;
[csDis setFloatValue:cstep];
return self;
}
/* this too ... */
-togFit:sender
{
stepmode=[sender intValue];
return self;
}

/* test routine (not used) */
-rotest:sender
{
printf("1: %f\n",[self dihedral:5 :1 :3 :6]*DRC);
[self rotateAbout:.0873 :1 :3];
printf("2: %f\n",[self dihedral:5 :1 :3 :6]*DRC);
[self update:U_POS];
return self;
}

/* set amino acid dihedral angles */
-setOmega:(float)ang
{
int i;
float a,b;

a=ang*DRC;
for (i=0; i<mol[0].numaa; i++) {
	if (mol[0].am[i].cn==-1||(!mol[0].am[i].sel)) continue;
	mol[0].am[i].omega=a;
	b=[self dihedral:mol[0].am[i].c1 :mol[0].am[i].n :mol[0].am[i].cn
		:mol[0].am[i].cn2];
	[self rotateAbout:a-b :mol[0].am[i].n :mol[0].am[i].cn];	
}
[self centerMol:0];
return self;
}

-setPsi:(float)ang
{
int i;
float a,b;

a=ang+180.0;
if (a>180.0) a-=360.0;
a*=DRC;

for (i=0; i<mol[0].numaa; i++) {
	if (!mol[0].am[i].sel) continue;
	mol[0].am[i].psi=ang*DRC;
	b=[self dihedral:mol[0].am[i].n :mol[0].am[i].c1 :mol[0].am[i].c2
		:mol[0].am[i].o];
	[self rotateAbout:a-b :mol[0].am[i].c1 :mol[0].am[i].c2];	
}
[self centerMol:0];
return self;
}

-setPhi:(float)ang
{
int i;
float a,b;

a=ang+180.0;
if (a>180.0) a-=360.0;
a*=DRC;

for (i=0; i<mol[0].numaa; i++) {
	if (!mol[0].am[i].sel) continue;
	mol[0].am[i].phi=ang*DRC;
	b=[self dihedral:mol[0].am[i].c2 :mol[0].am[i].c1 :mol[0].am[i].n
		:mol[0].am[i].h];
	[self rotateAbout:a-b :mol[0].am[i].c1 :mol[0].am[i].n];	
}
[self centerMol:0];
return self;
}

/* take 4 atoms and calculate a dihedral angle. The math is right, but */
/* it would take too much space to explain it. */
-(float)dihedral:(int)n1 :(int)n2 :(int)n3 :(int)n4
{
float cb[3],ca[3],bd[3],bc[3],v1[3],v2[3],v3[3];
float t,t2;

sub(&mol[0].coord[n2][0],&mol[0].coord[n3][0],cb);
sub(&mol[0].coord[n3][0],&mol[0].coord[n2][0],bc);
sub(&mol[0].coord[n1][0],&mol[0].coord[n3][0],ca);
sub(&mol[0].coord[n4][0],&mol[0].coord[n2][0],bd);

cross(cb,ca,v1);
cross(bd,bc,v2);
cross(v1,v2,v3);
t=dot(v3,bc)/len(bc);
t2=dot(v1,v2);

return(atan2(t,t2));
}

/* calculate a bond angle */
-(float)bondAng:(int)n1 :(int)n2 :(int)n3;
{
float v1[3],v2[3],r;

sub(&mol[0].coord[n1][0],&mol[0].coord[n2][0],v1);
sub(&mol[0].coord[n3][0],&mol[0].coord[n2][0],v2);

r=dot(v1,v2)/(len(v1)*len(v2));
return (acos(r));
}

/* center the molecule on the screen */
- centerMol:(int)m
{
float x,y,z;
int i,j;

x=y=z=0.0;
for (i=0; i<mol[m].numa; i++) {
	x+=mol[m].coord[i][0];
	y+=mol[m].coord[i][1];
	z+=mol[m].coord[i][2];
}

x/=(float)mol[m].numa;
y/=(float)mol[m].numa;
z/=(float)mol[m].numa;

for (i=0; i<mol[m].numa; i++) {
	mol[m].coord[i][0]-=x;
	mol[m].coord[i][1]-=y;
	mol[m].coord[i][2]-=z;
}

i=mol[m].am[4].o;
j=mol[m].am[7].h;
x=sqrt(SQR(mol[m].coord[i][0]-mol[m].coord[j][0])+
	SQR(mol[m].coord[i][1]-mol[m].coord[j][1])+
	SQR(mol[m].coord[i][2]-mol[m].coord[j][2]));

i=mol[m].am[4].o;
j=mol[m].am[6].h;
y=sqrt(SQR(mol[m].coord[i][0]-mol[m].coord[j][0])+
	SQR(mol[m].coord[i][1]-mol[m].coord[j][1])+
	SQR(mol[m].coord[i][2]-mol[m].coord[j][2]));

[hba setFloatValue:x];
[hb3 setFloatValue:y];

[self update:U_POS];
return self;
}

/* set the view window to a variety of convenient sizes */
-set320:sender
{
[vWindow sizeWindow:324 :202];
return self;
}

-set640:sender
{
[vWindow sizeWindow:644 :402];
return self;
}

-set6404:sender
{
[vWindow sizeWindow:644 :482];
return self;
}

-set533:sender
{
[vWindow sizeWindow:537 :402];
return self;
}

/* do a dihedral rotation about the line connecting n1 and n2 */
-rotateAbout:(float)a :(int)n1 :(int)n2
{
Atom *atom;
int i;
RtPoint n,r;
float mx[9],m;

a= -a;
atom=mol[0].atom;
for (i=0; i<mol[0].numa; i++) atom[i].tag=0;
atom[n1].tag=1;
tagatoms(atom,n2);
atom[n1].tag=0;
atom[n2].tag=0;

n[0]=mol[0].coord[n2][0]-mol[0].coord[n1][0]; /* n==vec from n1 to n2 */
n[1]=mol[0].coord[n2][1]-mol[0].coord[n1][1];
n[2]=mol[0].coord[n2][2]-mol[0].coord[n1][2];
m=sqrt(SQR(n[0])+SQR(n[1])+SQR(n[2]));
n[0]/=m;	/* n is unit vector */
n[1]/=m;
n[2]/=m;

m=1.0-cos(a);	/* set up mx transformation matrix */
mx[0]=cos(a)+n[0]*n[0]*m;
mx[1]=n[0]*n[1]*m+n[2]*sin(a);
mx[2]=n[0]*n[2]*m-n[1]*sin(a);
mx[3]=n[0]*n[1]*m-n[2]*sin(a);
mx[4]=cos(a)+n[1]*n[1]*m;
mx[5]=n[2]*n[1]*m+n[0]*sin(a);
mx[6]=n[0]*n[2]*m+n[1]*sin(a);
mx[7]=n[1]*n[2]*m-n[0]*sin(a);
mx[8]=cos(a)+n[2]*n[2]*m;

for (i=0; i<mol[0].numa; i++) {
	if (atom[i].tag==0) continue;
	r[0]=mol[0].coord[i][0]-mol[0].coord[n1][0];
	r[1]=mol[0].coord[i][1]-mol[0].coord[n1][1];
	r[2]=mol[0].coord[i][2]-mol[0].coord[n1][2];

	mol[0].coord[i][0]=r[0]*mx[0]+r[1]*mx[1]+r[2]*mx[2]+mol[0].coord[n1][0];
	mol[0].coord[i][1]=r[0]*mx[3]+r[1]*mx[4]+r[2]*mx[5]+mol[0].coord[n1][1];
	mol[0].coord[i][2]=r[0]*mx[6]+r[1]*mx[7]+r[2]*mx[8]+mol[0].coord[n1][2];
	atom[i].tag=0;
}

return self;
}

/* select atoms from SelectView data Structure */
-setSelAtom:(struct SELDAT *)list
{
int i;

for (i=0; i<mol[0].numa; i++) {
	if (list[i].sel) mol[0].atom[i].sel=1;
	else mol[0].atom[i].sel=0;
}

[molView display];
return self;
}

/* recursive atom tagger for rotation routines */
void tagatoms(Atom *atom,int n)
{
int i,j;

for (i=0; i<atom[n].numb; i++) {
	j=atom[n].bnd[i];
	if (atom[j].tag==0) {
		atom[j].tag=1;
		tagatoms(atom,j);
	}
}
return;
}

/* dot product */
float dot(float *a,float *b)
{
return (a[0]*b[0]+a[1]*b[1]+a[2]*b[2]);
}

/* cross product */
void cross(float *a,float *b,float *c)
{
c[0]=a[1]*b[2]-a[2]*b[1];
c[1]=a[2]*b[0]-a[0]*b[2];
c[2]=a[0]*b[1]-a[1]*b[0];
return;
}

/* vector subtract */
void sub(float *a,float *b,float *c)
{
c[0]=a[0]-b[0];
c[1]=a[1]-b[1];
c[2]=a[2]-b[2];
return;
}

/* vector length */
float len(float *a)
{
return (sqrt(SQR(a[0])+SQR(a[1])+SQR(a[2])));
}

/* return a random float between lo and hi */
float frand(float lo, float hi)
{
    return (((float)random() / (float)MAXINT) * (hi - lo) + lo);
}

@end